Dynamical thermalization and vortex formation in stirred 2D Bose-Einstein condensates 
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We present a quantum mechanical treatment of the mechanical stirring of Bose-Einstein condensates using 
classical field techniques. In our approach the condensate and excited modes are described using a Hamiltonian 
classical field method in which the atom number and (rotating frame) energy are strictly conserved. We simulate 
a r = quasi-2D condensate perturbed by a rotating anisotropic trapping potential. Vacuum fluctuations in the 
initial state provide an irreducible mechanism for breaking the initial symmetries of the condensate and seeding 
the subsequent dynamical instability. Highly turbulent motion develops and we quantify the emergence of a 
rotating thermal component that provides the dissipation necessary for the nucleation and motional-damping of 
vortices in the condensate. Vortex lattice formation is not observed, rather the vortices assemble into a spatially 
disordered vortex liquid state. We discuss methods we have developed to identify the condensate in the presence 
of an irregular distribution of vortices, determine the thermodynamic parameters of the thermal component, and 
extract damping rates from the classical field trajectories. 



PACS numbers: 03.75.Kk, 03.75.Lm, 05.10.Gg, 47.32.C- 

I. INTRODUCTION 

The experimental observation of vortex lattices in Bose- 
Einstein condensates (BECs) |[T] |2l E g] g] 13 H has stimu- 
lated intense theoretical interest and debate over the mecha- 
nisms involved in vortex formation. While vortex lattices had 
previously been observed in superconductors and superfluid 
Helium, dilute gas BECs offer the possibility of a quantita- 
tive description of the formation process using tractable the- 
ory. The experiments can be divided broadly into two cat- 
egories: (i) A thermal gas rotating in an axially symmetric 
trap is evaporatively cooled and condenses into a rotating lat- 
tice |6|; (ii) A low temperature condensate (r <K Tc ) is rota- 
tionally stirred, and vortices then nucleate and form a vortex 
lattice im |2] [3] m |5] |2l • The first category of experiment is 
conceptually simplest, and the origin of dissipative processes 
is well understood i8ll9l [T0ll . The second category of experi- 
ment is the subject of this paper, and as observed by a number 
of groups, has more complicated phenomenology. We shall 
concentrate on a typical, representative example, in which 
the condensate is initially in the interacting ground state of 
a cylindrically symmetric harmonic trap at a temperature very 
close to r = 0. Stirring is implemented by distorting the trap 
to elliptical symmetry, and rotating it at the frequency of the 
condensate quadrupole mode. This excites a dynamical insta- 
bility in which the condensate is transformed into a turbulent 
state, which later evolves into a rotating state containing a reg- 
ular vortex lattice. 

This process of vortex lattice formation by stirring provides 
a unique testbed for dynamical theories of cold bosonic gases, 
as noted by other authors yjLJ. Beginning from a well charac- 
terised initial state, the application of a simply defined experi- 
mental procedure causes a violent transformation into a highly 
excited and disordered state, from which an ordered state sub- 
sequently emerges. The workhorse of BEC theory, the Gross- 
Pitaevskii equation (GPE), can provide a good dynamical de- 



scription of r = condensates, and has given a very good 
understanding of the initiation of the dynamical instability 
IfTTl [121 [T3l [T4ll . The subsequent evolution into a highly ex- 
cited state is clearly beyond the validity regime of pure GP 
theory. Tsubota et al. [15] recognised the need for a new the- 
oretical framework after noting that dissipation is required to 
evolve an initially non-rotating GPE ground state into a state 
in equilibrium with the stirrer They postulated that a rotating 
thermal cloud develops to provide this dissipation (see also 
Ref 1 8 1), and modelled its effects by adding a phenomenolog- 
ical damping term to the GPE |16|, later justified 111] on the 
basis of the formalism of Jackson and Zaremba |17|. Further 
calculations employing the GPE were performed by Parker et 
al. IfTSl [191, and these 2D simulations gave rise to crystal- 
lization of a low energy vortex lattice, in conflict with earlier 
simulations of the GPE in 2D |20, 21 1 that did not exhibit lat- 
tice formation. It is of interest to identify the origin of this 
discrepancy and to determine the nature of the final state to be 
expected in stirring experiments of 2D Bose gases. 

In this paper we apply classical field theory, a represen- 
tation of many-body physics capable of modelling the dy- 
namical behaviour of both the condensate and excited states 
of the matter field. The formalism for classical field the- 
ory has been developed recently in a number of papers 
|22l mm |25l Ell [23, and it has been used to describe 
several phenomena where thermal fluctuations play a critical 
role in determining the dynamical behaviour of the system, 
including: the condensation transition UP. .27 J . effects of crit- 
ical fluctuations on the transition temperature |[28l . activation 
of vortices in quasi-2D systems |[29l[30ll3TI . collapse dynam- 
ics of attractive Bose gases |[32l . and coherence properties of 
spinor condensates Il33ll . Lobo et al. Il34ll have used a ver- 
sion of classical field theory to describe vortex lattice forma- 
tion by stirring, and their results indicate that the thermal field 
component generated provides the damping necessary for the 
lattice to form. Here we will apply an implementation of 
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classical field theory which conserves particle number and ro- 
tating frame energy, and uses a basis choice and numerical 
method that has appropriate boundary conditions and is free 
from numerical artifacts. This allows us to provide a detailed 
and quantitative description of the thermalization of the con- 
densate, and the nucleation and motional damping of vortices 
which is responsible for lattice formation. We find that the 
system relaxes to a state containing a disordered distribution 
of vortices, in thermal and rotational equilibrium with a cloud 
of thermal atoms generated by the stirring process. 

II. SYSTEM AND EQUATIONS OF MOTION 

The theory of the trapped Bose gas can be written in terms 
of the second quantized Hamiltonian for the field operator 

n 

where the operators satisfy bosonic commutation relations 
[an,a]n\ - 6„,n, and the modes 0„(x) are any complete or- 
thonormal basis. In the s-wave interaction regime the Hamil- 
tonian for the system is given by 

// = J flf^x 4'^(x)i/,p4'(x) 

-Hy ^ t/^x >I'^(x)4'^(x)4'(x)4'(x) (2) 

where U - Anfp-ajm, a is the scattering length and m is the 
atomic mass. The single particle Hamiltonian is generally of 
the form H^p - //"p -i- ye(x) -i- 6V(jl, t) where 

< = + Vo(x) - nL„ (3) 

and we have expressed the single particle Hamiltonian in a 
frame rotating at angular frequency Q. about the z-axis and 
decomposed the external trapping potential into 

y(x, f) = yo(x) + y,(x) + svi^, t). (4) 

We have singled out Vq(x) in Eqs. ([3]l, (|4]) as it will be used 
below to define the single particle basis for the system. V^ix) 
is the remaining time independent potential, and any time de- 
pendent potentials present are represented by 6V(x, t). 

A. The system 

We consider a condensate initially at temperature T - Q 
and confined in a cylindrically symmetric harmonic trap with 
trapping potential 

yo(x) = (m/2) [c4 + /) + ^2^2} . (5) 

The trap is deformed at time f = into an ellipse in the xy- 
plane which rotates around the z-axis with angular frequency 



Q. In the frame rotating at frequency Q, 6V - and the 
trapping potential takes the form 

y(x,f) = yo(x) + y.(x), (6) 

where 

y,(x) = ancjjiy^ - x\ (7) 

In this paper we set the trap eccentricity to e s 0.025 
which is sufficient to excite the quadrupole instability at Q ~ 
(jJr/ V2 IIT2irT3l . In this work we will always drive the system 
within the hydrodynamic resonance, choosing Q. = O.lSaJr- 
The sudden turn-on of the rotating potential at fixed e, Q, 
means that the rotating frame energy (the quantum average 
of Eq. (j2]i) is a conserved quantity of the motion. 

For computational reasons we now restrict our attention to 
an effective two dimensional system. Such systems have been 
experimentally realised (e.g., ||35]| ) by setting the trapping fre- 
quency o), sufficiently high that no modes are excited in the 
z direction. A 2D description of our system is thus valid pro- 
vided that Eg] 

fi + kT ^ fiLo,, (8) 

where yu and T are the system chemical potential and tem- 
perature respectively. Provided that the oscillator length I. = 
^hlniLj, corresponding to the z confinement greatly exceeds 
the s-wave scattering length a, the scattering is well described 
by the usual 3D contact potential lf36l . and we obtain an effec- 
tive two-dimensional equation in which the collisional inter- 
action parameter is simply rescaled to f/ao — 2 ^Jlnfr' a I mU 
ifTOl [36l . In this work we consider ^^Na atoms confined in a 
strongly oblate trap, with trapping frequencies {u)r, <jl>z) - Itiy. 
(10, 2000)rad/s. The s-wave scattering length is a = 2.75nm 
and the oscillator length L - 468nm. Having chosen appropri- 
ate parameters, in the next section we make the dimensional 
reduction formally rigorous by introducing an energy cutoff 
into the theory which excludes all excited z-axis modes. 



B. Classical field tlieory 

In projected classical field theory a cutoff energy is in- 
troduced to separate the system into a low energy region (L) 
with single particle energies satisfying e„ < Er, and a high en- 
ergy region (H) containing the remaining modes. Following 
ifTOl l26l I37II we shall refer to these regions as the condensate 
and non-condensate bands, respectively. This separation is 
motivated by the expectation that low energy modes will be 
significantly occupied. Where possible it is most convenient 
to carry out the separation in the basis which diagonalizes the 
single particle Hamiltonian of the system since at high ener- 
gies (~ Er) the many-body Hamiltonian (Eq. |2]i is approxi- 
mately diagonal in this basis (i.e., at such an energy a cutoff 
imposed at a single particle energy is approximately parallel 
to one imposed at an energy of the interacting system). 

For certain potentials an appropriate numerical quadrature 
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method exists which also greatly expedites simulations. Any 
numerical representation of field theory on a computer is also 
subject to a restriction of the modes. In what follows we pro- 
vide an explicit connection between the cutoff in our theory 
and our numerical implementation. 

We begin by introducing the projected field operator 



= 'p4'(x), 



(9) 



where the projector is defined as 

rm = ^ <f>„(x) f d'y 0:(y)/(y). (10) 

neL ^ 

The summation is over all modes satisfying e„ < Eji_, where 
H^p(f>„(x) = e„(p„(x) defines the modes in terms of a time in- 
dependent single particle Hamiltonian in the rotating frame. 
Choosing Er < tiUy, all excited z-axis modes are excluded 
from this set. In the s-wave limit, the effective low energy 
Hamiltonian for the system is then given by 



eff 



J flf^x iA"l"(x)H,p,A(x) + ■^«/'"'"(x)0-"^(x)f (X)f (X) (11) 



We now obtain a classical theory by making the replacement 
0-(x) i//(x) in //eff to give 



HcF = j d^x iA*(x)HspiA(x) 



(12) 



This approach is based on the high occupation condition, 
where it is known that classical behaviour is recovered for 
highly occupied bosonic modes since commutators become 
relatively unimportant Il23ll . We then obtain the equation of 
motion for the classical field via projected functional differen- 
tiation of the classical field Hamiltonian 1371 : 



ifi 



diffix) ^ 6HcF 
dt 6i//*{x) ' 



(13) 



The resulting equation of motion, the projected Gross- 
Pitaevskii equation (PGPE) 



in"^ = r {(//sp + f/2Dl«A(x)l') <A(x)} , 



(14) 



satisfies several identities governing the evolution of total 
number N - J d^x \if/(x)\^, rotating frame energy Hqf, and 
angular momentum L- |38 1: 



dN 
dt 
dT, 
dt 

dHcF 
dt 



--L,V(x)+-lm 
n n 



j d^X r(x) (y,(x) + f/2D|(^(x)|') Q {L,i/,ix)] 



2Q 



-Im 



j d^X .A*(x)(y.(x) + f/2Dl-A(x)l')e{L,lA(x)) 



(15) 
(16) 
(17) 



where (3 = 1 - !P is the projector orthogonal to f, and the 
bar represents spatial averaging: A - J d^x ij/*{x)At]/{x). The 
terms involving the imaginary part of an integral are boundary 
terms arising in the projected classical field theory. Although 
we work in the frame where 5V{x, t) = 0, Eqs. ( 14 1-( 17 1 hold 



for any partitioning of the time independent potential into ad- 
ditive parts V{){x) and V^ix) and therefore the precise choice 
of Vo(x) defining our projector in Eq 



and dL^jdt - {-i/fi)L.V^(x) as required. 

We emphasize that for certain potentials there exists an 
exact numerical quadrature method for evolving the PGPE 
which implements the projection operator to machine preci- 
sion for a given energy cutoff. Such a method exists for our 
choice of V'o(x) 1 10] and consequently there is an exact numer- 
ical correspondence between the formal properties of the clas- 



10 1 is in principle arbi- sical field theory expressed in Eqs. ( 15 i-( 17 1 and the results of 



trary. However, in order that we recover the formal properties 
of the continuous classical field theory, the boundary terms in 
Eqs. ( [T6| ), ([TtJi must be made to vanish. This is ensured by 
expanding il/{x) in eigenstates of L,, which is equivalent to re- 
quiring the projection operator to have rotational symmetry: 
[^,1,] = [Q,L,] = 0. We then have (3{L-(A(x)) = and we 
see that our choice of Voix) = mciyjir^/l given in Sec. II A 
generates the appropriate rotating frame classical field theory. 
The remaining potential V^{x) = ema>j(y^ - x^)/2 is not diag- 
onal in our representation, but since it is only a weak pertur- 
bation the many -body Hamiltonian will remain approximately 
diagonal in the eigenstates of //^p near the cutoff energy Er. 
With this choice of basis, we recover dN/dt = dHcF/dt = 0, 



PGPE simulation. PGPE evolution via numerical quadrature 
in the basis of L. eigenstates is thus free of boundary term 
artifacts, and we have obtained an appropriate Hamiltonian 
classical field theory in the rotating frame. 



C. Projected truncated Wigner method 

We have thus far confined our discussion to pure classical 
field theory without recourse to phase-space methods. Having 
established the appropriately conserving formalism we also 
include vacuum noise in our initial conditions, following the 
prescription of the truncated Wigner (TW) method [i22il . The 



4 



classical field i/r is formally related to the field operator by 
the Wigner function phase-space representation ll22l l24l l26l 
[37l [39], and thus even at zero temperature the initial state of 
i/r in any trajectory contains a representation of the vacuum 
fluctuations in the form of classical (complex) Gaussian noise 
of mean population equal to 5 -quantum per mode. It is of 
concern, then, to choose appropriate basis modes to which the 
noise is added (not to be confused with the basis modes in 
which the system is propagated), so as to faithfully represent 
the ground state of the many-body system [22 1. 

While higher order approaches which take account of the 
mutual interaction between the condensate and Bogoliubov 
modes exist |40|, in this work our primary aim is to include 
the possibility of spontaneous processes in the initial condi- 
tion. It therefore suffices to populate the Bogoliubov modes 
orthogonal to the condensate mode with noise |[22l [411 l42l . 
The initial state is constructed from the GP ground state of the 
symmetric trap in the lab frame, 0o(x), as 

i/r(x) = 0o(x) + Yj «/x)«; + Vj(x)a*j, (18) 
j 

where (uj, vj) are the Bogoliubov modes orthogonal to <;*o(x), 
Uj - sjrij + l/2(^j + iT]j)l V2, and for our system the thermal 
population nj is zero. The independent Gaussian distributed 
variables satisfy -rfj - rji^j = and ^i^j = jfiffj = 6ij. 

Population of the quasiparticle basis constructed in this 
manner ensures that all surface modes of the condensate (that 
are resolvable within the condensate band) are effectively 
seeded by noise. The dynamically unstable excitations of the 
condensate in the presence of the rotating trap anisotropy are 
among those seeded, and the noise introduced here thus plays 
a role entirely analogous to that played by vacuum electro- 
magnetic field fluctuations in triggering the spontaneous de- 
cay of a two-level atom |43 |. All symmetries of the mean- 
field state are thus broken |44| before the exponential growth 
of non-condensed field density B31 occurring during the dy- 
namical instability. By contrast an exact GP ground-state of 
the isotropic trap possesses circular rotational symmetry, and 
thus the exact evolution of the state under the GPE with an el- 
liptic potential would retain unbroken twofold rotational sym- 
metry for all time. 

In any practical calculation, with finite-precision arith- 
metic, numerical error eventually leads to the breaking of this 
symmetry, and researchers propagating the GPE from such an 
initial state have supplemented their numerics with additional 
procedures to speed this process |i8] [18] [191 . In this work we 
consider an irreducible source of symmetry breaking: vacuum 
fluctuations. 

There are also more subtle reasons for using a projec- 
tive method and for our choice of basis: (i) Preservation of 
symmetries. — It is natural in modelling the field theory of par- 
ticles interacting in free space, to consider the problem on a 
discrete spatial lattice B6l . effectively imposing a momentum 
cutoff on the system. However, the inclusion of an external 
confining potential breaks the translational symmetry of the 
system, and momentum is no longer a well-conserved quan- 
tity appropriate for defining a cutoff. A projective method al- 



lows us to implement a finite dimensional field theory which 
best respects the remaining symmetries of the confined sys- 
tem (e.g., the conservation of total field energy), (ii) Phase 
space discretization. — In modelling a homogeneous system 
on a discrete lattice, one expects to regain the full, contin- 
uum field theory as the lattice spacing tends to zero (i.e., as 
the momentum cutoff is raised to infinity). However, in a 
confined system, a natural discretization of the field theory 
is imposed by the quantization of energy levels in the con- 
fining potential. A cutoff defined in energy imposes both 
short-wavelength (ultraviolet) and long-wavelength (infrared) 
cutoffs consistently, and moreover, ensures that the number 
of modes spanning the intervening phase-space is optimally 
close to the true number present in the interacting system. Fur- 
thermore, the choice of an energy cutoff allows us to choose 
a phase-space appropriate to the rotating frame-the frame in 
which final equilibrium is achieved. This has immediate con- 
sequences for the distribution of thermal energy in the sys- 
tem. A non-optimal discretization, such as a Cartesian grid 
or planewave model of a trapped system, leads to somewhat 
arbitrary results for the number of modes, and consequently 
the temperature of the final equilibrium state. As we will 
see, our equilibrium states still have a cutoff dependent tem- 
perature, but the numerical prefactor arising from our choice 
of basis is close to optimal. A complete resolution of this 
problem requires a more general theory that includes above- 
cutoff corrections [47 1 . (iii) Boundary conditions. — Methods 
commonly used such as Fourier spectral methods ll34ll . and 
the Crank-Nicholson method \4E] can suffer from pathologies 
arising from periodic boundary conditions, wave vector alias- 
ing [49 [, or discretization error in the calculation of deriva- 
tives; both methods are non-projective and as such do not im- 
pose a consistent energy cutoff. For example, a model of the 
homogeneous Bose gas implemented using the Fourier spec- 
tral method will always suffer aliasing issues for any modes 
with momentum greater than half the maximum representable 
on the grid [50|. These issues are resolved by the use of a 
projective method. The system boundary is defined by the 
formally energy-restricted theory, and the numerically propa- 
gated basis modes and their coupling are unambiguously de- 
fined in direct correspondence to this theory. 



III. NUMERICAL IMPLEMENTATION 

We introduce now the dimensionless units we use in our im- 
plementation of the PGPE. We use harmonic oscillator units 
{{?, CO, E, t}) related to SI units ({r, a>, E, f )) by the expressions 
r - rro, u) - cjcor, E - ETlu),-, t - toj^^ where the radial 
oscillator length ro - yjTilmcOr. In reporting our results we 
shall often refer to times in units of the trap cycle (abbrevi- 
ated eye), i.e. the period of the radial oscillator potential, 
Tosc = Injci),-. The energy cutoff E^ separating the conden- 
sate and non-condensate bands is conveniently expressed in 
terms of yUi, and is typically chosen to be Zsr ~ 2 - 3/ijUi. This 
choice is made to be high enough so that: (i) there are suf- 
ficient available states to give a reasonable representation of 
thermalized atoms, and; (ii) the excited states at the cutoff are 



5 



only weakly affected by the mean-field interaction. However, 
ZiR must not be too large: the energy scale of our low en- 
ergy effective field theory must be well below that associated 
with the effective range of the interatomic potential, so that 
the contact s-wave potential remains valid [37J. Furthermore, 
the number of modes in our theory must be kept as low as pos- 
sible, so that the added vacuum noise does not dominate the 
real population |25 1. We examine the sensitivity of our system 
to the choice of cutoff in Sec. lIXI 



A. PGPE method 

Here we briefly discuss our implementation of the PGPE 
in terms of rotating frame harmonic oscillator states. A more 
thorough discussion of the method is presented in [lOJ . In our 



dimensionless units the PGPE (Eq. ( [T4| ) becomes 
.dip 



V2 ^2 



i^^'P{{-^ + '^{l-eco^29) + iQ.de + A\^\']ip}, (19) 
at 2 2 

where the dimensionless effective interaction strength A - 
U2Dlr\fujL>r- We proceed as in fTOl by expanding 41 as 



iA(x,f) = ^c„,(f)?„,(r,0), 
I",') 



(20) 



anisotropy 



HM)^-- d0 
^ Jo 



2n 



/-*0C 

Jo 



(r, ff)f- cos(20)iA(r, 9). (26) 



We implement the calculation of this term by means of a com- 
bined Gauss-Laguerre-Fourier quadrature rule similar to that 
introduced in 1 10] to calculate the F„i{iff), and refer the reader 
to Appendix | VIII| for the details. 



The PGPE as written in Eq. (24i immediately lends itself 
to numerical integration using the interaction picture method 
II5TI |52 I with respect to the isotropic single particle Hamil- 
tonian //'J,. Furthermore the absence of noise terms allows 

of the method. This 



us to use the adaptive-RK variant 
affords us explicit control over the level of truncation error 
arising during the integration. In practice we choose the ac- 
curacy of our integrator to be such that the relative change in 
field normalisation is < 4 x 10"^ per step taken by the integra- 
tor, and the simulations presented here all have total fractional 
normalisation changes < 2 x 10""* over their duration (~ 10^ 
trap cycles). Most importantly, the change in rotating frame 
energy is of magnitude commensurate with that of that change 
in normalisation, so that the truncation error represents a loss 
of population from the system as a whole, rather than a prefer- 
ential removal of high-energy components as in ifTSll . and the 
relaxation of the condensate band field is thus due to internal 
damping processes only. 



over the Laguerre-Gaussian modes 



(21) 



which diagonalize the isotropic single particle (non- 
interacting) Hamiltonian //'p 



i/°p?„,(r, 0) = [-^ + ^ + indg]Y„,(r, 6) = £„",?„,(?, 6), (22) 



with eigenvalues = 2n + \l\ - Q.I + \. The energy cutoff 
defined by the projector V requires the expansion of Eq. (20 1 
to exclude all terms except those in oscillator modes ?„i(r, 9) 
for which < E^, i.e. those modes satisfying 

2n + \l\-Q.l+ \ < £r. (23) 

This yields the equation of motion 

= E%, + AFM) + H„,W. (24) 
The projection of the GPE nonlinearity 



F„,(0-)= do 
Jo 



/-*0O 

-rd-rrj-r, 

Jo 



0m-r,ef;f,(-r,0), (25) 



is evaluated as in ifTOll . and this equation of motion (Eq. (24 1) 
differs only from the deterministic portion of the SGPE pre- 
sented there by the inclusion of the projection of the trap 



IV. SIMULATION RESULTS 

We present here the results of a simulation with initial 
chemical potential /ii = lAHo),-, whose response to the stir- 
ring is representative of systems throughout the range 4fi(L),- < 



jii < 2QhLij,-. Using the criteria discussed in Sec. Ill we set the 
condensate band cutoff Eji_ - 3/ii = 4-2fiLo,-. The upper limit 
of the range investigated (yUi ~ 2QTuOr) is set by computational 
limitations: at the fixed rotation frequency Q. - Q.lStOr the 
size of the Gauss-Laguerre (GL) basis scales as A4 ~ -E^, and 
the corresponding computational load scales as 0{E^ ifTOl . so 
that simulations rapidly become numerically expensive. For 
the cutoff /sr - A2h<jJr the GL basis consists of M - 2028 
modes. At the lower end of the range (yui < lOfiw^) the simula- 
tion results differ somewhat both in the thermalization process 
leading to vortex nucleation, and the behaviour of the vortex 
array, due to the reduced mean-field effects and small numbers 
of vortices respectively. These differences will be discussed in 
Sec. lVTAl 

We will focus primarily on the case of /ii = I4ha>r, and 
note for comparison that a pure condensate (i.e., ground GP 
eigenstate) with the same number of particles rotating at Q = 
0.75w,- would contain a lattice of 18 vortices. The response 
to the stirring is illustrated in the sequence of density distri- 
bution plots shown in Fig. [T] and exhibits the following key 
features. 

Dynamical instability and formation of thermal cloud. — 
The initial state (Fig. [TJ^)) quickly becomes elongated, with 
its long axis oscillating irregularly relative to the major axis 
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FIG. 1: (Color online) Classical field density at representative times during the evolution. Shown are (a) initial condition: lab. frame ground 
state plus vacuum noise, (b) ejection of material during the dynamical instability, (c) nucleation of vortices at the interface between the 
condensate and the non-condensate material, (d) penetration of vortices into the condensate bulk, (e) non-equilibrium state with rapid vortex 
motion, and (f) equilibrium state. Vortices are indicated by and are shown only where the surrounding density of the fluid exceeds some 
threshold value. Parameters of the trajectory are given in the text. 



of the trap. The quadrupole oscillations are dynamically un- 
stable to the stirring perturbation lfT2l [T3l IT4l . and since all 
the Bogoliubov modes of the condensate are effectively popu- 
lated due to the representation of quantum noise in our simu- 
lation, some of the fluctuations grow exponentially. The effect 
is clearly seen in Fig. [ijb), where matter streams off the tips 
of the elliptically deformed condensate in the direction of the 
trap rotation. Material is then ejected more or less continu- 
ously, forming a ring about the condensate, with the latter os- 
cillating in its motion and shape. The ring soon becomes tur- 
bulent and diffuse (Fig. [TJc)), losing coherence and forming 
a thermal cloud whose characterisation (Sec. |V A) is a central 
theme of this paper. In Sec. |V B 2] we show from an analysis of 
the angular momentum and density distribution that by f = 50 
eye. the outer part of the cloud (beyond r » 7ro) rotates as a 
normal fluid in rotational equilibrium with the drive. 

Vortex nucleation. — After the thermal cloud has formed, 
surface oscillations of the condensate grow in magnitude 
(Fig. JTTd)), in accordance with the prediction of Williams et 
ai, ||9] that such oscillations are unstable in the presence of 
a thermal cloud. These fluctuations form a transition region 
between the central condensate bulk, characterised by com- 
paratively smooth density and phase distributions (see also 
Fig. [Sja)), and an outer region of thermally occupied highly 
energetic excitations. Long-lived vortices are nucleated in this 
transition region, and then begin to penetrate into the edges 
of the condensate region. Initially the vortices penetrate only 



small distances r « ro with rapid and irregular motion, and 
visibly lag the rotation of the drive. Penckwitt et al. 18] have 
previously given a description of this vortex nucleation and 
penetration using a model in which a rotating thermal cloud 
was included phenomenologically. 

Vortex array formation. — The nucleated vortices remain 
near the edge of the condensate for a considerable time be- 
fore they begin to penetrate substantially into the condensate. 
Eventually vortices pass through the central region (by f ^ 75 
eye. in this simulation), with trajectories that are largely in- 
dependent, except when two vortices closely approach each 
other ll54ll . The density of vortices within the central region 
gradually increases (Fig.[T|e)), and the vortices begin to form 
a rather disorderly assembly, which still lags behind the drive 
in its overall rotation about the trap axis. The erratic motion 
of individual vortices and their differential rotation with the 
respect to drive then gradually slow, and they become more 
localised, until by f » 150 eye. they have formed into an array 
rotating at the same speed as the drive. This formation pro- 
cess can be interpreted as the damping of vortex motion by the 
mutual friction between them and the rotating thermal cloud 
155 , 56 1 . Alternatively, we can interpret the motion of the vor- 
tices as the result of low energy excitations of the underlying 
vortex lattice state f55\, which are gradually damped by their 
interaction with high energy excitations |57 , 58 , 59 1. The cou- 
pling of the high-energy modes to the low-energy excitations 
of the lattice state, drives their distribution towards equilib- 
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rium (W, '61] . The high energy modes comprising the ther- 
mal cloud equilibrate to a distribution with identifiable effec- 



tive thermodynamic parameters (see Sec. V Ai on a relatively 
short timescale (~ 100 eye). While the eff'ective temperature 
remains approximately fixed, the eff'ective chemical potential 
subsequently increases over a longer timescale (~ 3000 eye), 
as the distribution is shaped by its interaction with the low en- 
ergy excitations (including lattice excitations). By f a; 1000 
eye. the vortex motion reaches its equilibrium level, and the 
equilibrium populations of lattice excitations at this tempera- 
ture are such that the vortex distribution does not 'crystallize' 
into a rigid lattice. The damping of vortex motion is analysed 



quantitatively in Sec. V D 



V. ANALYSIS 

In this section we present the set of measurements we use 
to analyse the properties of the thermal cloud, and to iden- 
tify the remaining condensate fraction. We begin by making 
a simple estimation of the energy and temperature the thermal 
cloud will acquire during the lattice formation. We obtain the 
chemical potential of the thermal cloud by using a self con- 
sistent fitting procedure for the atomic distribution function, 
which also provides a more accurate measure of the temper- 
ature. We characterise the rotational properties of diff'erent 
spatial regions of the system, which provides evidence that 
during the vortex nucleation stage, the central part behaves as 
a superfluid, while the outer cloud behaves as a classical gas, 
rotating like a rigid body. 

The issue of condensate identification is a very important 
one, and is most commonly done in terms of spatial correla- 
tion functions, using ensemble averages of quantum mechan- 
ical operators Il62ll63]| . However, this is not well suited to the 
case of vortex lattice formation, due to varying vortex configu- 
rations within the ensemble, with correspondingly destructive 
effect on the spatial order of the system. In the classical field 
description of a finite temperature Bose gas, time averages are 
often substituted for ensemble averages, making use of the er- 
godic hypothesis |27|. We show, in Sec. VC that temporal 



correlation functions can be used to characterise the local co- 
herence of the field, and consequently determine the extent of 
the condensate mode and its eigenvalue. 

Following the first appearance of an identifiable vortex ar- 
ray, there is a long period over which it slowly relaxes to an 



equilibrium state. In Sec. V D we measure the motional damp- 
ing of the vortices and interpret it in terms of their interaction 
with the thermal component of the field |64|. Note that all 
results and figures in this section (Figs. 2]|9 1 correspond to the 
simulation shown in Fig.[T] 



A. Thermodynamic parameters 

1. Analytic predictions 

Since our system conserves both normalisation and energy 
in the rotating frame, some simple analytic predictions can 



be made about the development of the thermal component of 
the field, using the Thomas-Fermi (TF) approximation. In a 
frame rotating at angular velocity Q, the ground state is a vor- 
tex lattice which, neglecting the effect of the vortex cores, has 
Thomas-Fermi wavefunction 




where //q is the chemical potential of the rotating frame so- 
lution. From Eq. ( 27 1 we obtain the number of particles A^.^p 
and energy Ejp for the lattice state 



- 

JVjp - 



(28) 
(29) 



The corresponding quantities for the vortex free inertial 
frame ground state are given by setting Q = 0. We note that 
A'o = A'n=o has the same value in both the stationary and rotat- 
ing frames (since this solution is non-rotational). If we now 
compare a lattice state and a vortex free state, both fully con- 



densed and with equal occupation, we see from Eq. ( 28 1 that 
the chemical potentials are related by 



yUn = yUO -^1 - 02/w2^ 



(30) 



and hence from Eq. (29 1 that their rotating frame energies are 
related by 



TF — TF 



(31) 



The excess energy of the vortex free ground state over the 
lattice state of the same number of atoms. 



^^^Ej — ^^TF 



(32) 



can be significant, and for example is approximately for 
an angular velocity of f2 = 0.75w, . Thus in our stirring sce- 
nario, the rotating equilibrium state reached at the end of the 
process must have less atoms and less energy than the initial 
vortex free state, and the excess energy and atoms constitute a 
thermal cloud ll65l . The exact result for the excess energy of 
the classical field, obtained from the simulations, will depend 
upon the depletion of the condensate mode required to form 
the thermal cloud and the mutual interaction of the cloud with 
the condensate. For small thermal fractions this energy could 
be found by means of a calculation similar to that described 
in Ref. |44|. In the simplest approximation we assume the 
limit of a small thermal fraction, in which the energies of the 



condensate and thermal field are additive, with Eq. ( 32 1 ap 



proximating the excess energy the thermal field contains. We 
further assume that in equilibrium this energy will be classi- 
cally equipartitioned over weakly interacting harmonic oscil- 
lators, in the spirit of the Bogoliubov approximation |34|. The 
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equilibrium temperature can thus be predicted by 



in ESI 



M-l' 



(33) 



with Ai the condensate band multiplicity. For a given initial 
chemical potential we see therefore that the equilibrium tem- 
perature will be strongly dependent on the basis size. How- 
ever, we shall see in Section|Vl] Eq. ( [SO] ) provides a very good 
estimate of the final chemical potential reached by the field, 
which is essentially independent of cutoff. 



F(x,k) = 



1 



exp[(e(x,k)-Ai)/r]-l 



The semiclassical energy is 



e(x,K) 



2m 



+ V«/(x), 



(36) 



(37) 



with K the wave vector adjusted to the rotating frame fTOl, 
and the Hartree-Fock eff'ective potential in the rotating frame 
is 



2. Self-consistent fitting 

We expect a procedure such as the self-consistent Hartree- 
Fock approximation described in |66| to provide a good esti- 
mate for both fi and T. However in order to avoid performing 
the iterative procedure presented in |66| we use a simpler ap- 
proximation to this description, similar to that presented in 
^T\, but with some necessary modifications. The TF approx- 
imation for the condensate mode employed in |67 1 is a poor 
choice for the situation considered here, where the distribu- 
tion of vortices is irregular and the density of vortices is low. 
Secondly a straightforward application of the model of |67 | to 
our system leads to divergences, as the semiclassical density 
distribution of non-condensed atoms becomes 



nth(x) 



1 



gde 



-(V„r(x)+2£/2D(nc(>t)+«,h(x))-//)/A:7- 



(34) 



with yeif(x) - m(iol - Q?)r^/2 the centrifugally dilated trap- 
ping potential in the rotating frame. The Bose function gi{z) 
diverges logarithmically |68| as z — > 1, and thus the function 
nth(x) diverges where yeif(x) + 2f/2D(nc(x) -i- nth(x)) — > fi- Ne- 
glecting the mutual mean-field repulsion 2J/2D«th(x) of non- 
condensed atoms (as in |67|) would therefore leave us with 
divergences in nth(x)- We are thus led to develop a fitting 
function for the non-condensed density nth(x) only, and which 
therefore depends only on the total field density n(x), which 
we measure from the classical field trajectory. The inclusion 
of the mean-field interaction between non-condensed atoms 
ensures that our function remains regular so long as the (fitted) 
chemical potential is less than the minimum value of the effec- 
tive potential experienced by non-condensed atoms, a restric- 
tion which is not expected to constrain our fit in an unphysi- 
cal manner In this way we avoid both having to distinguish 
between condensed and non-condensed material in forming 
the mean-field potential, and having to iterate a self-consistent 
model to convergence. Our approach to determining // and T 
is as follows: The one-body Wigner function corresponding 
to a Bose field 0(x) is defined 



F(x,k) = 



dyifix + |)0(x ■ 



(35) 



and we assume here the non-condensate atoms in the con- 
densate band to be well described by such a distribution, for 
which we further assume the approximate semiclassical form 



yf/(x) = miojf - n'Y- + 2f/2D«(x), 



(38) 



with the total field density the sum of the condensate and non- 
condensate (thermal) contributions «(x) = nc(x) + nth(x). The 
factor of two here is the strength of the interaction between 
pairs of atoms when at least one of the pair is non-condensed 
[66, 67 69|. The classical field limit of Eq. (36i is obtained 



when F(x, k) » 1, which occurs when the parameter (e(x, k)- 
fj.)/kT « 1, and is given by its first order expansion in this 
parameter. 



Fc(x,k) 



kT 



e(x, k) - fi 



(39) 



It is important to note that although Eq. p9\ is obtained here 
as the high occupation approximation to the bosonic distribu- 
tion Eq. ( [36] l, (describing the classical equipartitioning of en- 
ergy over phase-space cells of volume l//i^), we expect it to 
apply to the equilibrium state of a weakly interacting classical 
field system even if the level occupations are < 1 |24. 34, 71 1 . 
With the local rotating frame wave vector cutoff defined by 
n^K^(xf/2m = maxl^R - yf/(x),0) HOl, we find for the 
spatial distribution 



n^^{x,f,,T) 



-I 



'^"f'" K dK 



In 



27T 



Fcix,K) 



(40) 



< rtp, with rtp 



^j2E^lm{ul - Vf-) the semiclassi- 
cal turning point of the condensate band in the rotating frame, 
and /IdB = ^j2T:tP- 1 mk^T the thermal de Broglie wavelength. 
This expression constitutes an appropriate fitting function for 
the non-condensed component of the classical field, which 
requires no explicit knowledge of the distribution of con- 
densed atoms ll72l . To obtain an estimate of and T for 
our system, we sample the azimuthally averaged field den- 
sity «(r) = Jj^ do n(r, 6) at A^, equally spaced times over 
some period f e [fo - T/2,to + T/2], and average over the 
sample times to obtain (n(r)) = Yjfli n{r,ti)INf From this we 
construct the effective potential 



<y^ff(r)) = m(w2 _ ^2yii2 + 2t/2D<«('-)). 



(41) 
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FIG. 2: (Color online) Fitting procedure for the thermodynamic pa- 
rameters. Shown are: the (time-averaged) radial density distribu- 
tion times radius (solid line), the effective potential experienced by 
non-condensate atoms times radius (dashed line), and the fitted dis- 
tribution of non-condensed atoms (dash-dot line). The dotted line 
indicates the classical turning point corresponding to energy cutolf 
E^. The data shown corresponds to the period t = 3000 - 3010 eye. 
(Fig.Qf)). 
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FIG. 3: (Color online) Evolution of the effective (a) temperature and 
(b) chemical potential of the non-condensed atoms, as measured by 
the semiclassical fitting procedure. 



B. Rotational parameters 



We then take as our fitting function 



1 



«f>;//,r)= — In 



(42) 



As our fitting function represents only the non-condensed 
component our field, we are obliged to avoid including sig- 
nificant condensate density in our fit to («(r)). We therefore 
perform our fit of njjj(r) to («(r)) over the domain r € [r_, rd, 
with r_ the location of the occurrence of the minimum value of 
the effective potential, i.e., {V™(r-)) = min{(y^^(r))), which 
should approximately mark the peak of the non-condensed 
fraction and thus the boundary of the condensate [67 ). Within 
the central region r < r_, the fitted density n^(r) thus decays 



monotonically as r 
Fig.|2] 



0. An example of such a fit is shown in 



After the initial strongly non-equilibrium dynamics follow- 
ing the dynamical instability (i.e. by f a; 400 eye), the aver- 
aged field densities («(r)) are fit well by the function given in 



Eq. (42 1, and we conclude that the higher energy components 
of the field have reached a quasistatic equilibrium. We fol- 
low thereafter the evolution of the temperature and chemical 
potential with time, which we present in Fig. |3] The tem- 
perature (Fig. ^a)) appears to be essentially constant within 
the accuracy of the measurement procedure, while the chem- 
ical potential (Fig. |3jb)) shows a definite upward trend over 
the course of the field's evolution. In Sec. VC2 we compare 



the chemical potential of the cloud determined from our fit- 
ting procedure, with the chemical potential of the condensate, 
which we extract by considering the temporal frequency com- 
ponents present in the field. 



A characteristic feature of superfluids is that they resist any 
attempt to impart a rotation to them, and they acquire angu- 
lar momentum through the nucleation of vortices |73], which 
endow the superfluid with quantized flow circulation (see e.g. 
Il56l ). Thus while a classical fluid in equilibrium with a con- 
tainer rotating at angular velocity Q. satisfies the equality 



where the classical moment of inertia is defined 



(43) 



(44) 



a superfluid in equilibrium with such a container will in gen- 
eral fail to do so. That is to say, the moment of inertia per 
particle 



1 = 



(45) 



will not equal the classical value I^/N in such a state. We note 
that for pure superfluids, such as a dilute Bose gas at T = 
(which is well described as a ground state of the GP equa- 
tion), the steady state angular momentum is typically much 
less than the classical value determined by its mass distribu- 
tion ll74ll . As the rotation frequency Q of the GP state in- 
creases, the (areal) density of vortices increases [75 , 76l, and 
is given to leading order by the Feynman relation [56] 



Pv = 



otQ 



(46) 



with corrections arising due to the lattice inhomogeneity and 
the discrete nature of the vortex array llTTll . As the rotation 
frequency and thus the vortex density increase, both the clas- 
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sical and quantum moments of inertia increase also, and the 
difference between the two moments vanishes in the limit of 
rapid rotations O — > ll76l . 



1. Vortex density 

To monitor the increase in vortex density during the sys- 
tem evolution, we count the vortices of positive rotation sense 
occurring within a circular counting region of radius equal to 
the TF radius of the initial state, /?xf - '<j2jdi/mcijj.. Due to 
the rotational dilation of the condensate in the rotating frame 
this region lies within the central bulk of the condensate, and 
in this manner we attempt to avoid including in the count the 
short-lived vortices constantly created and annihilated at the 
condensate boundary. At equally spaced times over a period 
of 4 trap cycles about time we obtain the average number 
of vortices in the counting radius, {n^{ti)), and from this we 
determine the average vortex density in the counting region 



Pv(f/) 



<«v(f,)) 



(47) 



^^v^^;.-,>;.,%.^._-v'^v■-''^■^>■•'•■'•■'■■''^^^^ 
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In Fi g.|4| a) we plot Eq. ( 47 1 as a fraction of the prediction of 
Eq. ( |46| l. During the dynamical instability and initial nucle- 
ation of vortices the vortex densities measured in this manner 
are spuriously high, due to the counting of phase defects in 
the turbulent non-equilibrium fluid rather than long-lived vor- 
tices in the condensate bulk. After this initial period, i.e. from 
t X! 150 eye. onwards, we see a gradual increase in pv as the 
condensate slowly relaxes and admits more vortices into its 
interior, and approaches rotational equilibrium with the trap. 
By f a: 6000 eye, the density appears to have essentially satu- 
rated at a level pv ~ O.SSp^. Such a value seems reasonable for 
the rotation rate and condensate size considered here, which 
are such that inhomogeneity and discreteness effects will be 
significant ffl\ . Moreover, it seems plausible that the high 
degree of lattice excitation here increases the energetic cost of 
vortices above that assumed in the mean-field description of 

Gil. 



2. Local rotational properties 

In order to characterise the rotation of the outer cloud 
formed from ejected material as well as that of the central 
bulk of the field we measure localised expectation values, i.e. 
expectation values over a restricted spatial region. In this way 
we can compare the mass distribution and angular momen- 
tum of particular regions, in order to characterise the localised 
properties of the field. We define the expectation value of an 
operator O on spatial domain D by 



{O). 



Jv 



(48) 



where the O will be either region ^ or S illustrated in Fig. |5] 
Focussing our attention on the outer annulus S, we calculate 



FIG. 4: (Color online) Rotational response of the field. Shown are: 
(a) Vortex density, with dash-dot lines indicating Pv/Pv = 0.7 and 
Pv/Pv = 1 for reference, (b) classical (dashed line) and quantum 
moments of inertia, and (c) ratio of quantum and classical moments 
of inertia. Insets show the behaviour of the same quantities over 
a longer timescale. Note that the initial L. is finite due to the ini- 
tial vacuum occupation. For reference dotted lines in (a-c) separate 
phases of the evolution discussed in the text. 



the classical moment of inertia of material in the annulus 



(/c)s = m{r 



(49) 



and the angular momentum of this material {L^}s. From mea- 
surements over the period f = 50 - 51 eye. we find that 



Q. 



(50) 



to within 4%, indicating that the cloud in this region is rotat- 
ing as a normal fluid in rotational equilibrium with the drive. 
Averaging samples of these expectation values over a period 
of a trap cycle, we find for the time averages {{L,}s/Q.} - 
1.025{(/c)s). 

Turning our attention to the central disc we find that 
over the same period, {(L,)g)/Q = 0.080{(/c)s). The cen- 
tral condensate bulk at this time thus remains approximately 
stable against vortex nucleation, but possesses some small an- 
gular momentum due to its shape oscillations, in contrast to 
the ejected material which has lost its superfluid character and 
has come to rotational equilibrium with the drive. 
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broken description as the coherent fraction of the field 
(Acoh = <iA), 



(51) 



10' 



hi 



FIG. 5: (Color online) Regions used for the evaluation of angular 
momenta and moments of inertia. Disc lies entirely within the 
condensate, where the flow is (initially) irrotational. Annulus S (be- 
tween the solid black line and the dashed line indicating the classical 
turning point of the condensate band) contains only non-condensed, 
normal fluid. The field density shown is that at r = 50 eye. 



3. Global rotational properties 

We consider now the rotational properties of the entire field, 
by evaluating expectation values as above over the full xy- 
plane, (i.e. in Eq. (48 1 we set D — > R^). A measure of the 



field's rotation is given by the angular momentum L~. We plot 
the ratio ///c to characterise the response of the atomic field 
to the imposed rotation, as both the field's angular momentum 
and its mass distribution change with time. In Fig. |4fb) we 
plot the evolution of these quantities as time evolves and in 
Fig.Qc) we plot the evolution of their ratio. 

Comparing the evolution of these quantities and that of 
the field's density distribution (c.f. Fig.[T]), we identify three 
roughly distinct phases of the system evolution: (i) Excitation 
of unstable surface mode oscillations accompanied by perma- 
nent increase in angular momentum (c.f Fig.[TJb)), (ii) further 
increase in angular momentum and Ijl^ as vortices are nucle- 
ated into the condensate bulk (Fig. [Tfc-e)), and (iii) gradual 
approach of the field to rotational equilibrium with the drive 
(Fig.[TJf)). Note that although the finite temperature equilibria 
here contain a significant normal fluid component, the quan- 
tum moment of inertia is still suppressed below the classical 
value due to the presence of the superfluid component. 



C. Temporal analysis 

1. Issues of condensate identification 

A central concern in classical field simulations is the iden- 
tification of the condensate mode and its occupation. In situa- 
tions such as 1 78 79 , 80 1 where the collisional dynamics dom- 
inate, and dynamics of interest occur on a short timescale, the 
condensate mode can be identified within a t/(l) symmetry- 



where the expectation value in the many-body state (■ ■ ■ ) can 
be formally related to moments of the classical field ESIED- 

However, on longer timescales the classical field is ex- 
pected to evolve to a thermal distribution Il23l l25l . and so 
we abandon the formal operator correspondences and exam- 
ine the classical statistics of the trajectory ensemble directly. 
Furthermore, the chaotic nature of the system means that tra- 
jectories which are initially close in phase-space will, in gen- 
eral, diverge as it evolves. Thus even in an idealised situation 
in which distinct trajectories each contain a highly occupied 
mode undergoing coherent phase rotation, the phase of this 
mode at a particular time t will vary between trajectories, and 
the identification of the condensate mode with the field opera- 
tor mean yields a condensate fraction which decays spuriously 
as the neighbouring trajectories dephase. This is a natural 
consequence of the symmetry-broken condensate definition, 
and occurs even for small non-condensate fractions ll42l [82l . 
however the dephasing is intuitively expected to occur more 
rapidly in systems with a greater incoherent fraction. 

However, global gauge symmetry breaking has no effect on 
the one body density matrix 



p(x,x') = <f'"(x)iA(x')), 



(52) 



or its classical field analogue 1271 . In a strongly thermal 
system, the ergodic hypothesis may be employed to approxi- 
mate such an ensemble expectation value by a time average. 
Blakie and coworkers |27| have used the Penrose-Onsager 
||83l definition of the condensate mode (the single eigenmode 
of p(x,x') with an occupation comparable to the total (mean) 
number of particles), to extract a condensate fraction from in- 
dividual classical field trajectories, in equilibrium and quasi- 
equilibrium scenarios. 

However, the Penrose-Onsager criterion is not appropri- 
ate for the current situation, which exhibits strongly non- 
equilibrium behaviour. As noted in [63.1 , in such situations 
the Penrose-Onsager definition may fail to correctly describe 
the amount of condensation in the system. Further complex- 
ity is added by the presence of a non-trivial phase structure 
such as that of a condensate containing an array of vortices. 
Indeed even at T = in an axially symmetric trap, the ro- 
tating frame ground state is in general highly degenerate due 
to the spontaneous breaking of the rotational symmetry of the 
many-body Hamiltonian 184J. The corresponding one-body 
density matrix thus fails to describe a single highly occupied 
mode. While in our case the axial symmetry of the Hamilto- 
nian is explicitly broken by the potential anisotropy ye(x), dis- 
tinct configurations of vortices distinguish states of the system 
in a way which cannot be cancelled by global phase rotations. 
One-body density matrices calculated from either ensemble 
or ergodic time averages thus describe a non-simple or: frag- 
mented condensate state (see |63|). As this behaviour might 
be expected of the true one-body density matrix in such a situ- 
ation, we regard this to be a lack in generality of the Penrose- 
Onsager criterion itself. A quantification of the condensate 
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population or mode at a single time would therefore require a 
description in terms of higher order correlation functions. 



2. Frequency distribution 

Previous authors studying equilibrium systems have used 
known quasiparticle bases to extract frequency distributions 
from classical field simulations |85, 86|. In a general non- 
equilibrium situation, with significant non-condensate popu- 
lations, suitable modes can not be defined. We therefore use 
an alternative strategy, and extract information about the dis- 
tribution of modes in the trap by analysing the frequencies 
present at particular points in position space. We define the 
power spectrum of the classical field i/r at position x about 
time fo 



//(x,t^;fo) = |g,^{tA(x,f))l', 
where (5^{/(0) denotes the Fourier coefficient 



gI{/(0) 



■io+T/2 



-T/2 



me 



"dt. 



(53) 



(54) 



We choose a sampling period of 10 trap cycles so as to be long 
compared to the timescales characterising the phase evolution 
of the condensate (rp ~ fi/fj), while being short compared to 
the that of the relaxation of the field (tr ~ 1000 eye.) ifTOl . 
and choose the sampling interval so as to resolve all frequen- 
cies present due to the combined single particle and mean- 
field evolutions (/iWmax = + We average this result 
over azimuthal angle 6j, to form H(r, uj, fo). In order to gauge 
the relative strength of the various frequency components at 
different radii within the classical region, we normalise the 
resulting data so that the total power at each radius Vj is the 
same. 

In Fig. |6|a), we plot the result of such a calculation over 
the time interval f = 9900 - 9910 trap cycles. For com- 
parison we have overlaid the centrifugally dilated potential 
Vc^X) - m((L)j. - Q^)r^/2 that would be seen by a classical par- 
ticle in the rotating frame, the cutoff frequency wr = E^^lfi, 
and the semiclassical turning point /"tp. Several notable fea- 
tures are present. At small radii we see a prominent peak in 
the power spectrum at w a; IQcOr, which we interpret as the 
coherent phase rotation frequency of the underlying conden- 
sate mode, i.e., the condensate eigenvalue yUc/fi. In Fig. |6|b) 
where the power spectrum is plotted at 3 particular radii, the 
peak in the spectrum obtained close to the trap centre (light 
gray/cyan line) can be seen clearly. In addition we see that 
frequencies above and below this are present, and we iden- 
tify these as due to the thermal occupation of the quasipar- 
ticle 'particle' and 'hole' modes respectively. At this same 
radial value (r - 0.066 Iro) a secondary peak is visible at 
u) =s 52(jL>r = (jUc -H Eb.)/Ti. An examination of the data re- 
veals that a smaller peak is also present at w a; -32(x),. = 
(pic - Ep}/h. This is an artifact of the projected method which 
occurs as follows: The Ai Bogoliubov modes which diago- 



nalize the Bogoliubov Hamiltonian in the condensate band in 
the presence of a stationary condensate mode (GP eigensolu- 
tion) can be considered as variational approximations to the 
lowest-lying 'true' Bogoliubov modes (obtained as the cutoff 
oo). The highest energy modes in the restricted quasi- 
particle basis are poor representations of the true modes, and, 
to maintain orthogonality of the set, are spuriously localised 
to the trap centre, and therefore have spuriously high ener- 
gies £k + fic ~ Ej^ + He due to the large mean-field interaction 
there. Therefore, assuming this behaviour of excitations to 
hold in a finite temperature equilibrium of the classical field, 
a peak is expected to occur at these frequencies at the trap cen- 
tre, where the density of these modes is comparatively large. 
However the population contained in the peak observed here 
is ~ 0.5% of the total thermal population and as such we do 
not expect it to qualitatively affect the relaxation process [87] . 
The black (blue) line in Fig.[6]]b) indicates the power spectrum 
at a larger radius where the peak is slightly less prominent. 
At larger radii still (r > 9ro) the energy distribution returns 
to approximately that of the non-interacting gas. Frequencies 
persist into the classically excluded region as a result of the 
familiar evanescent decay of energy eigenmodes into the po- 
tential wall. The effect of the mean-field interaction on the 
energy spectrum at this large radii, as can be seen in the dark 
gray (red) line in Fig. |6|b), is much less dramatic (e.g., only 
frequencies u > jidti are present), and the population of fre- 
quencies is only significant up to a small increment above the 
cutoff energy, which we interpret as the mean-field shift of the 
highest energy single-particle modes. We measure this shift to 
be a 2.5%, indicating that the cutoff we have chosen is suffi- 
ciently high to contain all modes significantly modified by the 



condensate's presence (see Sec. 11 B i. 

In order to follow the relaxation of the condensate, we 
measure the frequency Wp of maximum occupation in the 
power spectra, for power spectra evaluated over various peri- 
ods. Away from equilibrium this frequency does not necessar- 
ily correspond to the condensate eigenvalue, as many strong 
collective excitations are present. Nevertheless it provides a 
measure of the energy of condensate atoms, and we interpret 
it as representing an effective condensate chemical potential 
//p = h(jj-p. In Fig.|7]we plot the condensate chemical poten- 
tial obtained in this way and also the chemical potential of the 
thermal cloud obtained using our fitting procedure (Section 
V A2| l for comparison. We see that the condensate chemical 
potential reduces with time, as interactions with thermally oc- 
cupied modes damp its excitations. By f a 4000 eye, the two 
chemical potentials appear to have reached diffusive equilib- 
rium, to the accuracy of our fitting procedure for yu and T . 



3. Temporal coherence and condensate identification 

The power spectrum discussed above shows clearly that 
a condensate is present within the classical field, presenting 
as its signature the presence of a single highly occupied fre- 
quency, broadened by its interactions with the other modes in 
the system. The presence of this narrow frequency spike is 
indicative of the temporal coherence of this condensate mode. 
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FIG. 6: (Color online) (a) Azimuthally averaged power spectral density (PSD). Dashed lines indicate the cutoff energy, corresponding classical 
turning point and (centrifugally dilated) trapping potential, (b) PSD traces at particular radii. The black (blue) and dark grey (red) lines 
correspond to radii r = 3.1894ro and r = 11.9575ro respectively. The light grey (cyan) line corresponds to the smallest measured radius, 
(c) Autocorrelation function (absolute value) obtained from the power spectrum. The dashed line indicates the classical turning point of the 
condensate band, (d) Autocorrelation magnitude at radii corresponding to Fig.|6|b). Data corresponds to the period t = 9900 - 9910 eye. 
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FIG. 7: (Color online) Effective chemical potential of the condensate 
and thermal cloud. The data displayed is obtained from 10 trap cycle 
long samples spaced at 100 trap cycle intervals. Inset: Data with 
samples spaced at 10 trap cycle intervals over the initial decay (rise) 

of ;Up (yUth). 



and we will use the local autocorrelation function 

G(x,T;to)^WlrHH(x,cj;to)}, (55) 

of the field at point x, at lag t relative to the time to to quantify 
the local temporal coherence of the field. To do so we employ 
the Wiener-Khinchin theorem B9l to evaluate the local auto- 
correlation function. 

In practice we take the discrete Fourier transform of 
Hirii, cj; to) to obtain the azimuthally averaged autocorrelation 
G(rk, t; to), and normalise it so that |G(r,-, t = 0, fo)| - 1 for all 
radii r,-. We will use the modulus | ■ ■ ■ | to characterise the tem- 
poral coherence of the classical field, as a function of position, 
and take the timescale for its decay to represent the timescale 
over which the field remains temporally coherent at radius r^^, 
near time to- Illustrative results are shown in Figs. |6|c) and 
(d). 

At the centre of the trap (light gray/cyan line in Fig.|6|d)), 
the peak of the autocoiTelation function is broad and decays 
monotonically with increasing |r|. At larger radii (r » 2ro - 
5ro), this peak possesses a nontrivial structure (black/blue 
line), with \G\ vanishing and then increasing again as t varies. 
Analyzing the trajectories of vortices during the sampling pe- 



riod we conclude that this structure is due to the passage of 
vortices through these radii. We identify the width of the cen- 
tral peak (Fig. |6|d)) at a radius r, as (twice) the correlation 
time of the classical field at this radius. At the trap centre, the 
correlation time ~ 2.2 eye, indicating that despite the pres- 
ence of highly energetic thermal excitations and mode mixing 
in the central region, the temporal coherence of the condensate 
mode is significant. We note also that a small ridge appears in 
\G(t)\ for T close to zero, and interpret this as representing the 
correlations of the thermal component of the field (see also 

HMD- 

The width of the broad envelope (Fig.|6|d)) decreases with 
increasing radius around r - Sro, where the strong peak in 
the power spectrum also reduces. Past this point the peak be- 
comes much narrower (dark grey/red line in Fig. |6|d)); the 
short correlation time of the outer cloud reflecting the broad 
range of energies characteristic of the thermally occupied ex- 
citations. Beyond the classical turning point, the evanescent 
nature of the basis modes yields a spurious amount of tempo- 
ral coherence and this manifests as an increasing correlation 
time past this point (a similar phenomenon was observed for 
spatial correlations in ||88l ). 

While the interpretation of the correlation time we calcu- 
late in this manner is complicated by the non-monotonic de- 
cay of \G(t)\ with increasing |t|, it nevertheless allows us to 
distinguish turbulent behaviour from thermal behaviour For 
example, in Fig.[8|a) and (b) we plot the phase profile of con- 
densate band at time t - 100 eye, and the corresponding cor- 
relation time respectively. The correlation time maintains a 
near-constant value of fc ~ 1.1 eye. as r increases until r x tq, 
where the first vortex (phase singularity in Fig. |8ja)) occurs. 
At this point the correlation time begins to decay steeply with 
increasing r. Throughout the region r < r^p the phase struc- 
ture is complicated, however we see that the correlation time 
of the central region is at least of order t^ ~ 0.5 trap cycles, 
indicating that despite the turbulent nature of the field here, 
interactions serve to maintain its coherence, and we identify 
the behaviour here as superfluid turbulence 1891 , in contrast 
to the outer region (r > 6ro), which has a very short correla- 
tion time and thus appears to be truly thermal material. This 
is in contrast to the interpretation presented by the authors of 
lfT8l[T9l , who considered the turbulence developed during the 
stirring process to be purely superfluidic (i.e., zero tempera- 
ture) in nature. 



D. Vortex motion 

As noted in Sec. |IV| the motion of vortices (as viewed in 
the rotating frame) is initially very rapid, and slows as time 
goes on. In order to quantify this motion and its slowing with 
time, we track vortex trajectories, monitoring the coordinates 
of all vortices within a circular region of radius equal to the 
Thomas-Fermi radius of the initial lab-frame condensate. Due 
to the rotational dilation of the vortex-containing condensate, 
this region always lies within the central bulk of the conden- 
sate, as defined by the temporal analysis described above. In 
general, vortices enter and leave this counting region as time 
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FIG. 8: (Color online) Plots of (a) the condensate band phase and 
(b) the corresponding correlation time as a function of radius, near 
time f = 100 eye. 



progresses. We are interested in the motion of vortices which 
persist in the counting region and can therefore be considered 
to exist in the condensate, rather than those which occur in the 
violently evolving condensate periphery. We therefore discard 
the trajectories of all vortices which do not remain in the re- 
gion for at least half the counting period of = 10 trap cycles 
about time f,. 

The motion of the remaining vortices is erratic on a length 
scale of order of the healing length, which is the manifestation 
of thermal core-filling in the classical field model. As this mo- 
tion occurs on a short timescale, we interpret it as the signature 
of high-energy excitations of the vortices. Our interest how- 
ever is in the low energy component of vortex motion, which 
undergoes damping as the field relaxes to equilibrium. We 
therefore apply a Gaussian frequency filter centred on w = 
to the trajectory components {x(t),y(t)} to remove their high- 
est temporal frequency components. We choose the width of 
this filter to be cr^ ~ 0.6ciJr, and so the frequency components 
preserved in the trajectories correspond to low energy excita- 
tions (cij ~ LOr). From the A^, coordinate pairs {xk,yk) of the 
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filtered vortex trajectory we extract the mean speed along the 
trajectory arc, 

1 r 1 '^''^ I 

Vj(ti) = — I lis ~ — V J(xi - ;c,„i)2 + (ji-yi-i)^. 

(56) 

We then average this quantity over the vortices included in the 
count, to find the mean speed per vortex v(f,) = y^, v,(f,)/A^v- 
In Fig. we plot this quantity as a function of time. We note 
firstly that the mean vortex speed decays dramatically during 
the initial measuring period f « 70 - 150 eye, as vortices 
enter the condensate with high velocities (with a substantial 
component opposing the trap rotation, as observed in the ro- 
tating frame), and their motion is damped heavily by their in- 
teraction with one another and with the thermal component. 
Subsequent to this, a secondary, slower damping is observed 
over the period t ^ 150 - 1100 eye, after which the vortex 
motion seems to have reached a lower limit set by thermal 
excitations. During this period the condensate contains nearly 
the number of vortices it contains in equilibrium, and we inter- 
pret the damping during this period as the damping of low en- 
ergy lattice excitations due to their interaction with thermally 
occupied, high energy modes. Taking the log of the speed 
data (Fig.|9]inset), we see that the decay during this period is 
approximately exponential, and performing a linear fit to the 
data over this period we extract a damping rate y ~ 6.3 x 10""* 
(eye.)-'. 

We note that the final distribution of vortices is disordered, 
with the vortices failing to become localised in a regular lat- 
tice. This result is consistent both with the considerations of 
Sec. VA and with pure GPE studies of stirring in 2D per- 
formed by Feder et al. [20 1, and Lundh et al. ll2Tll . The 
results presented here (and those of ||34l ) suggest that the tur- 
bulent states observed in Il20l 1211 are in fact finite tempera- 
ture ones beyond the validity of the GP model employed in 
those calculations. By contrast, the calculations of Parker et 
al. IfTSl [T9l produced fully crystallized vortex lattices in 2D, 
indicating that additional damping mechanisms were available 
to the system in their simulations. It is possible that these are 
attributable to the spatial diff'erencing technique 1 481 l90l used 
in those calculations, which is known to be less accurate (91) 
than the methods used in the present work, and Refs. Ii20ii21j . 



VI. DEPENDENCE ON SIMULATION PARAMETERS 

We compare now the behaviour of stirred condensates for 
different condensate sizes, i.e., for different values of the ini- 
tial chemical potential fi^. The field density and thus the non- 
linear effects of mean-field interaction increase as the chemi- 
cal potential is increased, and larger condensates in rotational 
equilibrium will also contain larger numbers of vortices. In 



Sec. VI A we assess the effect of varying initial chemical po- 
tential on our classical field solutions. 

It is also important in performing classical field calculations 
such as those presented here to consider the effect of the cutoff 





height on the results of simulations. In Appendix IX we com 



FIG. 9: (Color online) Mean vortex velocity measured from the clas- 
sical field trajectory, revealing the rapid slowing of vortices during 
the initial stage of nucleation (t x - 120 eye), and more grad- 
ual damping of the remaining collective excitations (t x 100 - 1000 
eye). Inset: Logarithm of the mean vortex velocity (scaled by the 
minimum velocity measured Vn,in)> and linear least-squares fit to the 
data over the damping period. 



pare results for simulations with fixed initial chemical poten- 
tial but varying cutoff heights, and quantify the effect of this 
purely technical parameter on the physical predictions of our 
model. 



A. Chemical potential 

The initial chemical potential impacts upon the system evo- 
lution both qualitatively, affecting the density profile of the 
condensate mode and thus the nature of the dynamical insta- 
bility leading to vortex nucleation, and quantitatively, deter- 
mining the equilibrium parameters of the stationary state of 
the classical field. The characteristics of the vortex array and 
dynamics of the vortices are also strongly dependent on the 
chemical potential. We consider these three aspects in turn. 



1. Effect on dynamical instability 

The chemical potential determines the strength of the mean- 
field repulsion influencing the condensate's shape. As the 
chemical potential becomes large, the condensate mode ap- 
proaches that described by the Thomas-Fermi approximation, 
whereas for smaller chemical potentials, the kinetic energy of 
the condensate becomes important, and the condensate bound- 
ary broadens. Consequently the visual distinction between the 
condensate mode and its incoherent excitations becomes less 
clear as the chemical potential is reduced. This affects the 
early evolution of the system and the onset of the dynami- 
cal instability in a pronounced manner For the largest chem- 
ical potentials considered (fi ~ 17 - IQTuOr), the collective 
(quadrupole) excitations and their breakdown are clearly vis- 
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FIG. 10: (Color online) Evolution of the field angular momentum 
for various initial chemical potentials. 



ible against a background of incoherent noise resulting from 
the initial vacuum occupation. For smaller chemical poten- 
tials, while quadrupole oscillations are still visible and the 
outer cloud becomes visibly more dense as time proceeds, in- 
dicating that material is indeed lost from the condensate, any 
ejection of material from the condensate is obscured by the 
blurred condensate boundary. Persistent 'ghost' vortices can 
be seen forming in the incoherent region about the conden- 
sate periphery very early in the evolution (f ^ 20 eye. for 
fii = AtiiOr). In Fig. 10 we plot the evolution of the field an- 



gular momentum per atom for three different initial chemical 
potentials. The most obvious feature is that larger condensates 
support more angular momentum per atom, due to increased 
mean-field effects. The evolution of the field angular momen- 
tum for the smallest condensate considered (//; - AHoj^) re- 
veals the onset of irreversible coupling of angular momentum 
into the field at f » 130 eye. Inspection of the field density 
evolution shows that vortices begin to be nucleated into the 
central region of the field at about this time. We therefore 
conclude that despite the absence of visible break-up of the 
condensate, the unstable quadrupole oscillations still serve to 
produce the thermal component required to allow vortices to 
enter the central condensate region ||92i . 



2. Thermodynamic parameters 
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FIG. 11: (Color online) Equilibrium temperatures and chemical po- 
tentials reached as a function of initial chemical potential, as deter- 
mined by the fitting procedure. Dotted and dash-dot lines indicate 
the predictions of Eq. l |30^ and Eq. \33\ respectively. Numbers indi- 
cate times at which the measurements were made (thousands of trap 
cycles) in cases where numerical results were constrained by time. 

' eye. 



Eq. (33 1. We note that despite its approximate (TF-limit) and 
asymptotic (zero thermal fraction) nature, the analytical esti- 
mate for fiis a reasonable prediction of the measured values. 
The measured value of yU is slightly higher than the TF pre- 
diction, and this discrepancy appears to increase with increas- 
ing while the measured temperature appears consistently 
smaller than the analytical prediction at small yUi, and exceeds 
it as the initial chemical potential is increased. However, the 
agreement seems very reasonable given the simplicity of the 



arguments presented in Section V A which neglect both the 
kinetic energy of vortices (TF approximation), the depletion 
of the condensate population and all other effects beyond a 
linear (Bogoliubov) description. We therefore conclude that 
our simple analysis captures the essential physics determin- 
ing the thermodynamic parameters of the equilibrium. That is 
to say, the temperature is determined by the necessity of the 
system to redistribute the excess energy of the initial state (as 
viewed in the rotating frame), and that no additional heating 
of the system is required for the system to migrate to a state 
containing vortices. 



As described in Section VA a simple analysis based on 



energetic considerations allows us to predict the equilibrium 
temperatures and chemical potentials of our classical field 
simulations, and the dependence of these parameters on the 
initial chemical potential and the condensate band multiplicity 
of our simulations. In Fig. 1 1 we plot these analytic predic- 



tions alongside the values obtained using the fitting procedure 
of Section lVAl The dotted line shows the estimated chemical 
potential of the (idealised) norm conserving, rotating frame 
condensate mode, jiQ - jUi -y/l - fi^/w^, while the dash-dot 
line indicates the estimated equilibrium temperature given by 



3. Vortex array structure 

Perhaps the most striking difference between vortex arrays 
in simulations of differing condensate populations is the num- 
ber of vortices present. As noted in Sec. VB the density 
of vortices is largely determined by the rotation rate of the 
condensate and is given to leading order by Eq. (46 1, and 



(within the TF approximation) the total number of vortices 
in an equilibrium condensate is approximately proportional to 
the chemical potential fx. 

The stability of a vortex array in equilibrium is determined 
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by the competition between hydrodynamic and energetic con- 
siderations |i56| which favour the crystalHzation of a rigid lat- 
tice, and the perturbing effect of thermal fluctuations in the 
atomic field. Close to zero temperature these fluctuations are 
interpreted as the thermal occupation of Bogoliubov excita- 
tions of the underlying vortex lattice state. As the temperature 
of a vortex lattice state is increased it eventually undergoes a 
transition to a disordered state 1 10| as the long-range order of 
the lattice is degraded by thermal fluctuations. 

After a long time period, f ~ 10"^ trap cycles, all our sim- 
ulated fields appear to describe states on the disordered side 
of this transition. All our solutions have vortices migrating 
within the condensate bulk, and leaving and entering through 
the condensate periphery. In the the smallest condensates, 
with the smallest vortex counts (i.e. // = (4, 5,6)/;^, with 
counts A^v ~ (4,5,6)), the vortex array is small enough that 
quasi-regular configurations of vortices may occur, in which 
the vortex positions appear to fluctuate about approximate 
equilibrium positions on a triangular (or square) lattice. Such 
configurations may persist for as long as ~ 10 trap cycles, but 
ultimately breakdown and give way to new configurations as 
vortices cycle in and out of the central, high density region 
of the field. We interpret this cycling behaviour as the classi- 
cal field 'sampling' different configurations during its ergodic 
evolution |93|. 

Larger vortex arrays appear to be prohibited from forming 
regular structures over scales larger than nearest-neighbour 
inter- vortex separations (the scale of the apparent order in the 
smaller condensates). As noted in Sec. V A the equilibrium 



(f) 



temperature attained by our simulated atomic fields increases 
linearly with the initial chemical potential, and so it seems 
unlikely that this behaviour would abate when the condensate 



size is increased further. Fig. 12 a-c) we plot the coordinate 
space densities of vortex arrays of differing sizes, correspond- 
ing to differing initial chemical potentials. Because of their 
small size compared to the extent of the condensate, vortex 
cores are difficult to experimentally image in situ, and so their 
presence is typically detected after free expansion of the atom 
cloud |94| by absorptive imaging techniques. For condensates 
in the TF regime the expansion is well described by a simple 
scaling of the position-space density (95^ "96] (though the vor- 
tex core radius grows in somewhat greater proportion ||94| ). 
The Beer-Lambert law ||97l for the absorption of an optical 
probe, and the assumption of a ground-state (Gaussian) pro- 
file in z yields for the transmitted intensity of the probe 



I(x,y) = ke 



(57) 



with /() the input probe intensity and y a constant which de- 
pends, in general, on the absorptance of the atomic medium. 
In Figs. [T2|d-f ) we plot simulated intensity profiles obtained 
using Eq. (|57[), where we have in each case simply set y - 
2/{|i/'|~)maxi with {|i/'P)niax the peak value of the density occur- 
ring in the distribution, so as to best represent the significant 
features. This process accentuates both the disparity in density 
between the condensate and the outer thermal cloud, and the 
density fluctuations in the central condensate region, in con- 
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FIG. 12: (Color online) (a-c) Coordinate space densities of the clas- 
sical fields with initial chemical potentials yU; = 5, 10, ISSWr at times 
t = 9900,9900,8800 eye. respectively, (d-f) Simulated absorption 
images generated from the same densities. 



them. The images show the thermally-distorted vortex array 
structure that might be measured in an experiment. The fact 
that the formation of rigid vortex lattices here is inhibited by 
thermal fluctuations, while ordered lattices have been formed 
experimentally at comparatively high temperatures fT,'4^, sug- 
gests that this inhibition of lattice order may be a manifes- 
tation of the well-known susceptibility of two-dimensional 
system to long-wavelength fluctuations |98|. Indeed (to our 
knowledge), no investigations of the effect of rotation on the 
BEC and BKT phases of Bose gases in quasi-2D trapping ge- 
ometries have been performed to date ||99]| . 



4. Damping of vortex motion 



While the vortex motion as measured by the mean veloc- 



trast to the logarithmic plots (Figs. 12 a-c)) which suppress 



ity defined in Sec. V D exhibits the same overall damping be- 
haviour in all our simulations, the damping appears to be gen- 
erally non-exponential and quite sensitive to the initial noise 
conditions. We have observed additional features (e.g. tran- 
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sient and oscillatory behaviours superposed with the damping) 
peculiar to each trajectory. It is therefore difficult to unam- 
biguously quantify the rates of motional damping in order to 
compare the dependence on initial chemical potential, with- 
out generating ensemble data for each parameter set, which 
would be a very heavy task numerically. We therefore plot in 
Fig. 13 a) only the behaviour of three particular trajectories. 



which indicate at a qualitative level that the vortex motion in 
the smaller condensates damps to equilibrium more quickly. 
We note however that the level of vortex motion that each 
trajectory damps to shows a clear dependence on the initial 
chemical potential. In Fig.[T3|b) we plot the mean velocity of 
vortices over the period t - 1900 - 2000 eye, at which point 
the motion in each simulation appears to have essentially re- 
laxed to its equilibrium level. The mean velocity so measured 
clearly increases with chemical potential. After comparing 
the results for varying cutoff at fixed initial chemical potential 
(not shown) we conclude that this is likely due to the increase 
in equilibrium temperature with increasing initial condensate 
size. 



VII. CONCLUSIONS 

We have carried out the first strictly Hamiltonian simula- 
tions of vortex nucleation in stirred Bose-Einstein conden- 
sates. Our approach is free from grid method artifacts such 
as aliasing and spurious damping at high momenta, enabling 
a controlled study of the thermalization and vortex nucleation 
within an explicitly Hamiltonian classical field theory. 

Vacuum symmetry breaking. — The importance of symmetry 
breaking has been highlighted in previous works |8, 18, 19]. 
Sampling of initial vacuum noise provides an irreducible 
mechanism for breaking the twofold rotational symmetry of 
the classical field solution. 

Thermalization. — Resonant excitation of the quadrupole 
instability evolves the system from a zero temperature ini- 
tial state through a dynamical thermalization phase in which a 
rotating thermal cloud forms. Subsequent vortex nucleation 
in the remaining condensate shows that the cloud provides 
the requisite dissipation to evolve the condensate subsystem 
toward a new quasi-equilibrium state in the rotating frame. 
This picture is consistent both with energetic constraints of a 
time independent rotating frame Hamiltonian and with previ- 
ous treatments based on introducing dissipation into the GP 
equation IfSlfTSl. 

Identifying the thermal cloud. — We have quantified the 
heating caused by dynamical instability through a range of 
measures. Spectral analysis of ergodic classical field evolu- 
tion has been used to identify the condensate chemical po- 
tential, and thermal cloud properties were extracted using a 
semiclassical fit at high energies. While the chemical poten- 
tials of the condensate and emerging thermal cloud are slower 
to equilibrate, the moments of inertia, temperature, and angu- 
lar momentum of the classical field are rapidly transformed by 
the instability to those of a rotating, heated Bose gas. 

Frustrated crystallization. — Vortices are seen to nucleate 
and enter the bulk of the condensate but a regular Abrikosov 
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FIG. 13: (Color online) (a) Vortex motional damping for 3 repre- 
sentative choices of initial chemical potential /jj. (b) Mean vortex 
motion (averaged over the period t = 1900 - 2000 eye. as a function 
of initial chemical potential /Ji. 



lattice does not form. Instead the vortices distribute through- 
out the condensate in spatially disordered vortex liquid state, 
consistent with substantial thermal excitation of an underly- 
ing regular vortex lattice state. In this respect our results are 
consistent with previous 2D GPE treatments of similar stirring 
configurations. 



We conclude that vortex nucleation by stirring is a finite 
temperature effect, which arises from dynamical thermaliza- 
tion when initiated from a zero temperature BEC. While we 
observe vortex nucleation in 2D, the asymptotic absence of 
lattice rigidity suggests that the temperature attained in ac- 
commodating the excess (rotating frame) energy of the initial 
state is such that vortex lattice crystallization is inhibited in 
this minimal, conserving stirring scenario in 2D. 
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where the are the roots of L^Jx) and the weights Wk are 
known constants 111 001 . We can thus write Eq. ([58]l as 



W, N„-l 

X cos{4n jl/Ng)i//{ V^,6ij), 



(62) 



where the effective weights - rnvk^k- Once the modes are 
precomputed at the Gauss points 



Pi = 1>„|,|(X,), 



(63) 



VIII. CALCULATION OF THE TRAP ANISOTROPY 



In this appendix we outline our numerical method for 
calculating the dimensionless perturbing potential V^{r) — 
-|r-cos(20). To calculate this term we follow IfTOll in ex- 
ploiting the properties of the Laguerre-Gaussian basis modes 
to evaluate this perturbation exactly using a Gauss-Laguerre 



quadrature rule. Changing variables Eq. (26 1 becomes 



e 

"4^ 



■Ik 



dOe'"" cos(20) 



(58) 



X I dxxi>„i(x)if/(y/x,6)7T, 
Jo 



where ^ni(x) = !'„/( V^, 0)e The integral to evaluate is 
therefore of form 



1 r^" 

47r Jo 



Jo 



dOe-"" cos(26) dxe-''Q(x,0), (59) 



where Q{x,0) is a polynomial in x and e'" of order deter- 
mined by the cutoff. The order of e'" in i/'(x) is constrained 
by Eq. (jgs} to be -L = -L(0) < I < l+iQ) = /+, where 
l±(n) - [{N - 2«)/(l + Q)], with [...] denoting the floor func- 
tion. Choosing the 6 grid so as to accurately compute all 

integrals ^ ^ d9e~"'^, where q is an integer in the range 

-(/_ + 1+ + 2) < q < (~L + 1+ + 2), the integral can be 
carried out exactly for the condensate band by discretizing 
6 as 6j = jAO = jln/Ng, for j = 0, 1, ■ ■ ■ , A^e - 1, where 
Ng = 1- +1+ + 3, to give 



■27T 



r. N„-l 

-'■«cos(20)=9^ J]< 



-ilnjllNu 



Nf, 



cos(26lj)2(x, Oj). 



i=0 



(60) 

The jc-integral is computed by noting that the polynomial part 
of the x-integrand when the angular integral is nonzero is of 
maximum order [A^/(l - Q.)] = l+, which can be computed 
exactly using a Gauss-Laguerre quadrature rule of order = 
[1+/2] + 1. The radial integral in Eq. (59 1 is 



dxQ(x,6)e-' ^y wkQ(xk,0), (61) 

ti 



the mixed quadrature rule for calculating the perturbation due 
to the trap anisotropy is: 



1 . Compute the radial transform 



lN/2] 

n=0 



(64) 



2. Construct the position field at the quadrature points '1'^:^ 
using the EFT, after zero padding xu to length Ne in the 
/ index, = pad,(^'«, A^e), 



N„-\ 



1=0 



3. Apply the perturbing potential 

Fkj = 'i'kjXkCos{29j). 

4. Compute the inverse FFT 



^ j=0 



-ilnjl/Ne 



(65) 



(66) 



(67) 



5. Calculate the Gauss-Laguerre quadrature 



(68) 



k=\ 



The orders of the quadrature rules used here are significantly 
smaller than those used for the evaluation of the nonlinear 
term (see ITOl ). and consequently the computational load in- 
crease due to the inclusion of the trap anisotropy over that of 
the base method is of order 25-50%, depending on the rotation 
rate and cutoff" height. 



IX. CUTOFF DEPENDENCE 

Of crucial importance in classical field simulations is the ef- 
fect of the basis size upon the system behaviour. As noted by 
other authors Il25ll34l . for simulations in which the developed 
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FIG. 14: (Color online) Field temperatures and chemical potentials 
reached as a function of cutoff. Dashed and dash-dot lines indicate 
the predictions of Eq. ^30\ and Eq. l[33| respectively. Measurement 
times are indicated as in Fig.|l 1| 




FIG. 15: (Color online) Decay of the effective condensate chemical 
potential (as defined in Sec. |V C| l for trajectories with initial chemical 
potential fi, = lOTia),. and three different values of the energy cutoff 
Er. 



thermal fraction is small, such that a Bogoliubov-level de- 
scription of the equilibrium non-condensate distribution is ap- 



propriate, the temperature is expected to be inversely propor- 
tional to the modes available for thermalization. It might be 
expected |34 | that this will influence relaxation rates, which 
depend strongly on the thermal occupation of system excita- 
tions (see e.g. I.57J ). 



5. Thermodynamic parameters 



In Fig. 14 we plot the temperatures and chemical potentials 
attained in simulations with //; = 20 and cutoff heights 
varying over the range — > 3//i. We note first that while 
the temperature shows a clear dependence on Er, the chem- 
ical potential appears relatively insensitive to it. The system 
therefore reaches a chemical potential closest to that result- 
ing from the idealised transformation of the condensate to a 
lattice state without loss of atoms from the condensate mode. 
The measured temperature decreases with increasing cutoff as 
expected, but differs from the T oc 1/(A1 - 1) behaviour pre- 
dicted by for a classical field in the Bogoliubov limit. We 
attribute this to the constraint imposed by the low cutoff. The 
cutoff here serves to increase the energy of natural modes in 
the system, by increasing their overlap with the condensate 
(c.f. the discussion in Sec. VCj, thus lowering the equilib- 
rium temperature of those modes. Indeed for = 2//, the 
classical turning point rtp ~ l.VrxF, where rjF ~ 7.8ro is the 
TF radius of the lattice state, and so the corruption the quasi- 
paiticle mode shapes may be quite pronounced. 



6. Damping 



As noted in Sec. VD it is difficult to extract quantitative 
comparisons of relaxation rates from the damping of vortex 
motion. Here we consider the relaxation of another quantity: 
the effective condensate chemical potential (power spectrum 



peak) /ip introduced in Sec. V C In Fig. 15 we show the decay 
of this quantity with time for three values of the cutoff Er. 
We observe that the damping is reasonably insensitive to the 
cutoff. A possible explanation for this behaviour is that as 
the effect of lowering the cutoff is to decrease the number of 
thermally occupied modes in the system and consequently to 
increase the temperature of these modes, that these two effects 
serve to approximately cancel one another in determining the 
rates of thermal damping. 
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